Abstract
This notebook is a template for evaluating COVID-19 hospitalization forecast submissions from COVIDhub. After inputting a set of parameters (forecasters, COVID signals, etc), the template yields a comprehensive report on the predictions of COVID forecasters as well as their performance compared to the ground truth. The visualizations generated by the template offer an intuitive way to compare the accuracy of forecasters across all US states.
\[\\[.4in]\]
Every week, forecasters submit their hospitalization predictions to COVID-19 ForecastHub. In this report, we rely on an AWS bucket that contains the estimates of a handful of signals (e.g., COVID death, cases, hospitalization, etc). Furthermore, the AWS server stores an array of evaluation metrics of these forecasts (e.g., Absolute Error, Weighted Interval Score, and 80% Coverage). Alternatively, the data can be retrieved from the publicly accessible covidcast and covideval APIs.
s3a <- get_bucket("forecast-eval", region = "us-east-2")
s3b <- get_bucket("forecasting-team-data")
scores1 <- s3readRDS("score_cards_state_hospitalizations.rds", s3a,
region = "us-east-2")
scores1 <- subset(
scores1,
forecaster %in% union(params$forecasters, "COVIDhub-baseline"))
# The Hub always makes forecasts on Monday (dofw 2)
# We move all forecast_dates to the following Monday and shorten the ahead
wday_shift <- function(x, base_dofw = 2) (base_dofw - wday(x)) %% 7
scores1 <- scores1 %>%
mutate(ahead = ahead - wday_shift(forecast_date),
forecast_date = forecast_date + wday_shift(forecast_date))
# need to pass in the right forecaster here
scores <- s3readRDS(params$aws_scores, s3b) %>%
magrittr::extract2("reformatted.buggy.matched.scorecards") %>%
select(ahead, geo_value, forecaster, forecast_date, target_end_date,
actual, wis, ae, cov_80, value_50, data_source, signal,
incidence_period)
# Logan's names are brutally long...
our_pred_dates <- unique(scores$forecast_date)
forecast_dates <- our_pred_dates
aheads <- unique(scores$ahead)
geo_values <- unique(scores$geo_value)
veggies <- c("asparagus", "broccoli", "chard", "daikon",
"escarole", "fennel", "garlic", "horseradish",
"jicama", "kohlrabi", "lettuce", "mushroom",
"nori", "okra", "parsnip", "radish", "squash",
"tomato", "watercress", "baseline1", "yuca")
names_table <- tibble(forecaster = unique(scores$forecaster), veggies = veggies)
scores <- left_join(scores, names_table) %>%
select(-forecaster) %>%
rename(forecaster = veggies)
results <- scores1 %>%
select(ahead, geo_value, forecaster, forecast_date, target_end_date,
actual, wis, ae, cov_80, value_50, data_source, signal,
incidence_period) %>%
bind_rows(scores) %>%
filter(forecast_date %in% forecast_dates,
ahead %in% aheads,
geo_value %in% geo_values)
The target forecast dates are:
2020-11-16, 2020-11-23, 2020-11-30, 2020-12-07, 2020-12-14, 2020-12-21, 2020-12-28, 2021-01-04, 2021-01-11, 2021-01-18, 2021-01-25, 2021-02-01, 2021-02-15, 2021-02-22, 2021-03-01, 2021-03-08, 2021-03-15, 2021-03-22, 2021-03-29, 2021-04-05, 2021-04-12, 2021-04-19, 2021-05-03, 2021-05-10, 2021-05-17, 2021-05-24, 2021-05-31, 2021-06-07, 2021-06-14, 2021-06-21, 2021-06-28, 2021-07-05, 2021-07-12, 2021-07-19, 2021-07-26, 2021-08-02, 2021-08-09, 2021-08-16, 2021-08-23, 2021-08-30, 2021-09-06, 2021-09-13, 2021-09-20, 2021-09-27
The template will compile data of the following forecasters:
COVIDhub-ensemble, USC-SI_kJalpha, MOBS-GLEAM_COVID, Karlen-pypm, JHUAPL-SLPHospEns.
For this analysis, all of Logan’s forecasters have been renamed:
kableExtra::kbl(names_table) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))
| forecaster | veggies |
|---|---|
| ______________________none dacouno zacouno neimacpas | asparagus |
| ______________________none dacouno zacouno seimacpas | broccoli |
| ______________________none oacouno zacouno neimacpas | chard |
| ______________________none oacouno zacouno seimacpas | daikon |
| ______________________none oareigo zareifo neimacpas | escarole |
| ______________________none oareigo zareino neimacpas | fennel |
| ______________________none oareino zareino neimacpas | garlic |
| ______________________none zacougo dacouno neimacnah | horseradish |
| ______________________none zacougo dacouno neiwinnah | jicama |
| ______________________none zacouno dacouno neimacnah | kohlrabi |
| ______________________none zacouno dacouno neiwinnah | lettuce |
| __________________never-up oareigo zareino neimacpas | mushroom |
| _________________always-up oareigo zareino neimacpas | nori |
| ________________always-low oareigo zareino neimacpas | okra |
| _______________always-down oareigo zareino neimacpas | parsnip |
| _____________always-steady oareigo zareino neimacpas | radish |
| ___marginal_cheating_u/sdl oareigo zareino neimacpas | squash |
| __nqcomb_logistic1_u/s/d/l oareigo zareino neimacpas | tomato |
| _marginal_cheating_u/s/d/l oareigo zareino neimacpas | watercress |
| baseline1 | baseline1 |
| marginal_logistic1_u/s/d/l oareigo zareino neimacpas | yuca |
\[\\[.07in]\]
Mean <- function(x) mean(x, na.rm = TRUE)
GeoMean <- function(x, offset = 0) exp(Mean(log(x + offset)))
facets.label = str_glue("{aheads} days ahead")
names(facets.label) = aheads
subtitle = sprintf("Forecasts made over %s to %s",
format(min(forecast_dates), "%B %d, %Y"),
format(max(forecast_dates), "%B %d, %Y"))
plot_wis <-
plot_canonical(results,
x = "forecast_date",
y = "wis",
aggr = GeoMean,
grp_vars = c("forecaster","ahead"),
facet_rows = "ahead", dots = FALSE,
base_forecaster = "COVIDhub-baseline") +
labs(title = subtitle,
x = "Forecast Dates",
y = "Geometric Mean WIS") +
#geom_point(aes(forecast_date, round(wis, digits = 2), color)), alpha = 0.05) +
facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead=facets.label)) +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
scale_y_log10() +
geom_hline(yintercept = 1, size = 1.5) +
scale_color_viridis_d() +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_wis, tooltip="text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))
plot_wis <-
plot_canonical(results,
x = "forecast_date",
y = "wis",
aggr = Mean,
grp_vars = c("forecaster","ahead"),
facet_rows = "ahead", dots = FALSE,
base_forecaster = "COVIDhub-baseline") +
labs(title = subtitle,
x = "Forecast Dates",
y = "Geometric Mean WIS") +
#geom_point(aes(forecast_date, round(wis, digits = 2), color)), alpha = 0.05) +
facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead=facets.label)) +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
scale_y_log10() +
geom_hline(yintercept = 1, size = 1.5) +
scale_color_viridis_d() +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_wis, tooltip="text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))
plot_wis_a <-
plot_canonical(results,
x = "ahead",
y = "wis",
aggr = GeoMean,
grp_vars = c("forecaster"),
dots = TRUE,
base_forecaster = "COVIDhub-baseline") +
labs(title = subtitle,
x = "Days ahead",
y = "Geometric Mean WIS") +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
geom_hline(yintercept = 1, size = 1.5) +
scale_y_log10() +
scale_color_viridis_d() +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_wis_a, tooltip="text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))
plot_wis_a <-
plot_canonical(results,
x = "ahead",
y = "wis",
aggr = Mean,
grp_vars = c("forecaster"),
dots = TRUE,
base_forecaster = "COVIDhub-baseline") +
labs(title = subtitle,
x = "Days ahead",
y = "Geometric Mean WIS") +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
geom_hline(yintercept = 1, size = 1.5) +
scale_y_log10() +
scale_color_viridis_d() +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_wis_a, tooltip="text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))
plot_cov80 <-
plot_canonical(results,
x = "forecast_date",
y = "cov_80",
aggr = mean,
grp_vars = c("forecaster","ahead"),
facet_rows = "ahead",
dots = FALSE) +
labs(title = subtitle, x= "Forecast date", y = "Mean Coverage 80") +
facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead = facets.label)) +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
scale_color_viridis_d() +
geom_hline(yintercept = 0.8, size = 1.5) +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_cov80, tooltip="text", height=800, width=1000)
plot_cov80_a <-
plot_canonical(results,
x = "ahead",
y = "cov_80",
aggr = mean,
grp_vars = "forecaster",
dots = TRUE) +
labs(title = subtitle, x= "Days ahead", y = "Mean Coverage 80") +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
scale_color_viridis_d() +
geom_hline(yintercept = 0.8, size = 1.5) +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_cov80_a, tooltip="text", height=800, width=1000)
library(sf)
results_intersect <- intersect_averagers(
scores, c("forecaster"), c("forecast_date", "geo_value")) %>%
select(c("ahead", "geo_value", "forecaster","forecast_date", "data_source", "signal","target_end_date","incidence_period","actual","wis","ae","cov_80"))
kl <- function(q, p = .8) p*log(p/q) + (1-p)*log((1-p)/(1-q))
maps <- results_intersect %>%
group_by(geo_value, forecaster) %>%
summarise(wis = Mean(wis),
cov_80 = Mean(cov_80),
kl_80 = kl(cov_80)) %>%
left_join(animalia::state_population, by = "geo_value") %>%
mutate(wis = wis / population * 1e5) %>%
pivot_longer(c("wis", "cov_80", "kl_80"), names_to = "score") %>%
group_by(score) %>%
mutate(time_value = Sys.Date(),
max = max(value), min = min(value)) %>%
group_by(forecaster, .add = TRUE)
keys <- maps %>% group_keys()
maps <- maps %>% group_split()
levs <- levels(maps[[1]]$score)
# for county prediction, set geo_type = "county"
maps <- purrr::map(maps,
~as.covidcast_signal(
.x, signal = .x$score[1],
data_source = .x$forecaster[1],
geo_type = "state"))
maps <- purrr::map2(
maps, keys$score,
~plot(.x,
choro_col = scales::viridis_pal()(3),
range = switch(.y,
cov_80 = c(0,1),
wis = c(.x$min[1], .x$max[1]),
kl_80 = c(0, .x$max[1]))))
cowplot::plot_grid(plotlist = maps[keys$score == "wis"], ncol = 3)
cowplot::plot_grid(plotlist = maps[keys$score == "cov_80"], ncol = 3)
cowplot::plot_grid(plotlist = maps[keys$score == "kl_80"], ncol = 3)